RPM_RPMA_Allos_CellTag notebook

Modifying original code

Original file: Fig3_4_6_ExtFig5-6-10_RPM_RPMA_Allos_CellTag.R

Only CellTag analysis included in this notebook.

Related to:

• Fig 3c-h
• Extended Data Fig 6b,c
• Extended Data Fig 7f

suppressPackageStartupMessages({
    library(Seurat)
    library(SeuratObject)
    library(SummarizedExperiment)
    library(ggplot2)
    library(ggpubr)
    library(tidyr)
    library(patchwork)
    library(viridis)
})

SAVEFIGS <- FALSE

TBO_seurat<-readRDS("../data/05_2025_RPM_RPMA_TBO_CellTag_Seurat_wSigs_FA_dpt_final.rds")
TBO_seurat

An object of class Seurat 
110982 features across 26618 samples within 2 assays 
Active assay: RNA (55491 features, 0 variable features)
 2 layers present: counts, data
 1 other assay present: norm
 2 dimensional reductions calculated: umap, fa

clones <- readRDS("../data/05_2025_RPM_RPMA_TBOAllo_CellTagClones_Onlyclones.rds")

my_colors <- c(
  "#E41A1C", # strong red
  "#377EB8", # medium blue
  "#4DAF4A", # green
  "#984EA3", # purple
  "#FF7F00", # orange
  "#FFFF33", # yellow
  "#A65628", # brown
  "#e7298a", # pink
  "#666666", # grey
  "#66C2A5", # teal
  "#FC8D62", # salmon
  "#8DA0CB", # soft blue
  "#E78AC3", # soft pink (different from 8)
  "#A6D854", # light green (but yellowish tint, not green)
  "#FFD92F", # lemon yellow
  "#E5C494", # light brown
  "#B3B3B3", # light grey
  "#1B9E77", # deep teal
  "#D95F02", # dark orange
  "#7570B3", # strong purple
  "#66A61E"  # olive green (NOT same green as before)
)

pheno_col <- c("brown2","darkorchid4","dodgerblue","#66A61E","orange","turquoise4","turquoise")
names(pheno_col) <- c("NE","NE/Neuronal","Neuronal","ATOH1","Tuft","Triple-Neg","Basal")

Fig 3c

ForceAtlas embedding of RPM and RPMA CellTagged allografted tumor cells

DimPlot(TBO_seurat, group.by='leiden_scVI_1.2', cols=my_colors, reduction='fa', label=TRUE, label.size=4) & 
    NoAxes() + theme(legend.position="none") # suppress legend

DimPlot(TBO_seurat, group.by='Pheno', cols=pheno_col, reduction='fa', label=FALSE, label.size=6) & 
    NoAxes()

# Extract coordinates from force-directed layout (or any layout)
dat <- Embeddings(TBO_seurat, reduction = "fa") %>% 
  as.data.frame() %>%
  mutate(Cluster = TBO_seurat$leiden_scVI_1.2,
         Genotype = TBO_seurat$Genotype)

# Set up your color vector (named)
my_colors_named <- setNames(my_colors, levels(TBO_seurat$leiden_scVI_1.2))

# Plot each Genotype with background cells greyed
plots <- lapply(unique(dat$Genotype), function(g) {
  dat$highlight <- ifelse(dat$Genotype == g, "highlight", "background")
  dat$Color <- ifelse(dat$Genotype == g, my_colors_named[dat$Cluster], "lightgrey")

  ggplot(dat, aes(x = FA_1, y = FA_2)) +
    geom_point(aes(color = Color), size = 0.5) +
    scale_color_identity() +
    ggtitle(g) +
    theme_void() +
    theme(
      plot.title = element_text(hjust = 0.5),
      legend.position = "none",
      plot.background = element_rect(fill = "transparent", color = NA)
    )
})

# Combine into one figure
a <- wrap_plots(plots, ncol = 2)
a

if(SAVEFIGS) ggsave("FAbyGeno.png", plot= a, width=6, height=3, dpi=300, bg = "transparent")

Visualizing Celltag/clone data

From table above, added whether clones were robust or not
Robust defined as >5 cells per clone post-QC.
CellTag metadata: clone information can also be found in Supplementary Table 4 of Ireland et al, 2025

Idents(TBO_seurat) <- 'Robust'

table(TBO_seurat@meta.data$Robust, TBO_seurat@meta.data$UnID)

        
         RPMA_Allo RPM_Allo3 RPM_Allo4 RPM_Allo_New
  No          6659      5244      6523         2227
  Robust      3501      1474       337          653

table(TBO_seurat@meta.data$Robust, TBO_seurat@meta.data$GenoCT)

        
         RPMA_CTpreCre RPM_CTpostCre RPM_CTpreCre
  No              6659         11767         2227
  Robust          3501          1811          653

Identify robust clones

table(clones@meta.data$Robust, clones@meta.data$Genotype)

        
          RPM RPMA
  Robust 2464 3501

Assess clones by leiden cluster

Fig. 3d

First, just look at bar graph, no particular order #

Idents(clones) <- 'leiden_scVI_1.2'

x <- table(clones@meta.data$CellTag_Clone,Idents(clones))
proportions <- as.data.frame(100*prop.table(x, margin = 1))

# proportions$Cluster
colnames(proportions)<-c("Cluster", "Sample", "Frequency")

# Stacked
p <- ggplot(proportions, aes(fill=Sample, y=Frequency, x=Cluster)) + 
  geom_bar(position="stack", stat="identity")

p + scale_fill_manual(values=my_colors) + 
  theme_bw()+ theme(axis.text.y = element_text(size=20), 
                    axis.text.x=element_text(size=14), axis.title.x =element_text(size=14), 
                    axis.title.y = element_text(size=18), legend.text = element_text(size=12), 
                    legend.title = element_text(size=18))+rotate_x_text(size=7,angle = 90)

dat <- p$data

Transform data to perform hierarchical (agglomerative) clustering

# Pivot to wide format: Cluster × Sample
mat <- dat %>%
  pivot_wider(names_from = Sample, values_from = Frequency, values_fill = 0) %>%
  tibble::column_to_rownames("Cluster") %>%
  as.matrix()

mat_norm <- prop.table(mat, margin = 1)  # normalize each row to sum to 1

Hierarchical (agglomerative) clustering of clones with default parameters in pheatmap()

hm <- pheatmap::pheatmap(t(mat), cutree_rows = 1, cutree_cols = 8, cellwidth = 5, 
                         cellheight = 5, fontsize = 8,
                         cluster_rows=FALSE, border_color=NA, 
                         color = colorRampPalette(c("darkturquoise","black","red2"))(30))

clone_order <- colnames(t(mat))[hm$tree_col$order]

Fig 3d

Clustered clones by proportions of cells in Leiden clusters

# Now, plot clones by Leiden, ordered, final # (Fig. 3d)
Idents(clones)<-'leiden_scVI_1.2'

x <- table(clones@meta.data$CellTag_Clone,Idents(clones))
proportions <- as.data.frame(100*prop.table(x, margin = 1))

colnames(proportions) <- c("Cluster", "Sample", "Frequency")
proportions$Cluster <- factor(proportions$Cluster, clone_order)

# Stacked
p <- ggplot(proportions, aes(fill=Sample, y=Frequency, x=Cluster)) + 
    geom_bar(position="stack", stat="identity")

p + scale_fill_manual(values=my_colors) + 
  theme_bw() + theme(axis.text.y = element_text(size=20), 
                     axis.text.x=element_text(size=14), axis.title.x =element_text(size=14), 
                     axis.title.y = element_text(size=18), legend.text = element_text(size=12), 
                     legend.title = element_text(size=18))+rotate_x_text(size=7,angle = 90) + 
    labs(x = NULL, y = "% of clone")

Idents(clones) <- 'Clone_Dynamics'
table(clones@meta.data$Clone_Dynamics)

Pattern_1 Pattern_2 Pattern_3 Pattern_4 Pattern_5 Unknown_1 Unknown_2 
      553       675       693      2797      1032       167        48 

Function for plotting clone dynamics

clone_patterns <- c('Pattern_1', 'Pattern_2', 'Pattern_3', 'Pattern_4', 'Pattern_5', 'Unknown_1', 'Unknown_2')
pattern_colors <- c('orange', 'green2', 'red', 'royalblue2', 'purple', 'gray40', 'black')
names(pattern_colors) <- clone_patterns

plotCloneDyn <- function(clone_pattern, pattern_color, tit = "") {
    # Get cells of interest
    ss <- subset(clones, idents = clone_pattern)
    highlighted_cells <- rownames(ss@meta.data)
    
    # Plot all cells in grey, highlight the pattern in your chosen color
    DimPlot(clones, 
            reduction = "fa", 
            cells.highlight = highlighted_cells,
            sizes.highlight = 0.25,
            cols.highlight = pattern_color,
            cols = "grey90", 
            pt.size = 0.1) + 
    ggtitle(tit) & NoLegend() & NoAxes()
}

Fig 3e

ForceAtlas embedding showing representative clones of each pattern with cells colored by pattern

cdp <- (lapply(1:7, function(i) plotCloneDyn(clone_patterns[i], pattern_colors[i], tit=clone_patterns[i])))
p_combined <- ggarrange(plotlist = cdp, ncol = 4, nrow = 2)
p_combined

getCloneCellIDs <- function(id) {
    cellIDs <- rownames(subset(clones,idents=c(id))@meta.data)
}
clone_cell_ids <- lapply(setNames(clone_patterns, clone_patterns), getCloneCellIDs)

plotPhenoWithBackground <- function(clone_pattern, pheno_col, title = "", seurat_obj = clones) {
    # Check if 'fa' reduction exists
    if (!"fa" %in% names(seurat_obj@reductions)) {
        stop("The 'fa' reduction is not found in the Seurat object.")
    }

    # Extract FA coordinates and metadata
    coords <- as.data.frame(Embeddings(seurat_obj, reduction = "fa"))
    coords$cell <- rownames(coords)

    # Get metadata
    meta <- seurat_obj@meta.data[, c("Clone_Dynamics", "Pheno")]
    meta$cell <- rownames(meta)

    # Merge coordinates and metadata
    dat <- merge(coords, meta, by = "cell")

    # Flag highlighted cells
    dat$highlight <- dat$Clone_Dynamics == clone_pattern
    dat$Pheno <- droplevels(dat$Pheno)

    # Split background and foreground
    bg_dat <- dat[!dat$highlight, ]
    fg_dat <- dat[dat$highlight, ]

    if (nrow(fg_dat) == 0) {
        warning(paste("No cells found for Clone_Dynamics pattern:", clone_pattern))
        return(ggplot() + ggtitle(paste(clone_pattern, "(no cells)")))
    }

    # Plot
    p <- ggplot() +
        geom_point(data = bg_dat, aes(x = FA_1, y = FA_2), color = "gray90", size = 0.2) +
        geom_point(data = fg_dat, aes(x = FA_1, y = FA_2, color = Pheno), size = 1) +
        scale_color_manual(values = pheno_col, drop = FALSE) +
        ggtitle(title) +
        theme_void() +
        theme(
            legend.position = "none",
            panel.background = element_rect(fill = "transparent", color = NA),
            plot.background = element_rect(fill = "transparent", color = NA)
        )

    return(p)
}

ForceAtlas embedding showing cells colored by SCLC phenotype

DimPlot(clones, group.by='Pheno', cols=pheno_col, reduction='fa', label=FALSE, label.size=6, pt.size = 0.05) & 
    NoAxes()

Fig 3f

ForceAtlas embedding showing representative clones of each pattern with cells colored by SCLC phenotype

plots <- lapply(clone_patterns, function(pat) plotPhenoWithBackground(pat, pheno_col, title = pat))
p_combined <- ggarrange(plotlist = plots, ncol = 4, nrow = 2)
p_combined

Visualize individual clones

Ext Data Fig. 6b,c

plotPatternDynByClone <- function(clone_pattern) {
    Idents(clones)<-'Clone_Dynamics'
    p1 <- subset(clones,idents=c(clone_pattern))
    
    pattern_color <- pattern_colors[clone_pattern]

    test <- as.data.frame(p1$CellTag_Clone)
    test$Barcodes <- rownames(test)
    # Group by CellTag_Clone
    test <- test %>% group_by(p1$CellTag_Clone)
    test <- dplyr::group_split(test)
    
    n_clones <- length(test)
    # Plot in for loop all RPM clones in Pattern 1
    plot_lst <- vector("list", length = n_clones)
    for (i in seq(n_clones)) {
      g <- DimPlot(TBO_seurat, group.by="CellTag_Clone", reduction='fa', order=TRUE, 
                   cells.highlight=test[[i]]$Barcodes,sizes.highlight=2, 
                   cols.highlight=pattern_color) +
          ggtitle(paste0(test[[i]]$`p1$CellTag_Clone`[1])) + 
          theme(plot.title = element_text(size = 8,face = "plain")) & NoLegend() & 
          NoAxes()
      plot_lst[[i]] <- g
    }
    
    # Combine multiple plots for output, as desired
    return(cowplot::plot_grid(plotlist = plot_lst, ncol=5))
}

pattern_plot_list <- lapply(clone_patterns, plotPatternDynByClone)

n_plots_per_pattern <- sapply(pattern_plot_list, function(x) length(x[['layers']]))

## calculate height of figure based on number of plots per pattern
plt_height <- 2.667 * ((n_plots_per_pattern %/% 5) + 1)

Ext Data Fig 6b,c

pattern_plot_list[[1]]

pattern_plot_list[[2]]

pattern_plot_list[[3]]

pattern_plot_list[[4]]

pattern_plot_list[[5]]

pattern_plot_list[[6]]

pattern_plot_list[[7]]

Visualize only clones matching in vivo in FA projection

Ext Data Fig. 7f

Idents(clones) <- "CellTag_Clone"
# Subset just the clones that match the in vitro, pattern 1
invitromatch <- subset(clones, idents=c("RPM_Clone_14","RPM_Clone_2","RPM_Clone_33","RPM_Clone_36","RPM_Clone_6"))

p1_cells <- colnames(invitromatch)
p2_cells <- colnames(subset(clones, idents="RPM_Clone_23"))
p5_cells <- colnames(subset(clones, idents="RPM_Clone_13"))

DimPlot(TBO_seurat, group.by="CellTag_Clone", reduction='fa', 
        order=TRUE, cells.highlight=p1_cells, sizes.highlight=2, 
        cols.highlight=c("orange")) + ggtitle("") & NoLegend() & NoAxes()

DimPlot(TBO_seurat, group.by="CellTag_Clone", reduction='fa', 
        order=TRUE, cells.highlight=p2_cells, sizes.highlight=2, 
        cols.highlight=pattern_colors["Pattern_2"]) + ggtitle("") & NoLegend() & NoAxes()

DimPlot(TBO_seurat, group.by="CellTag_Clone", reduction='fa', 
        order=TRUE, cells.highlight=p5_cells, sizes.highlight=2, 
        cols.highlight=pattern_colors["Pattern_5"]) + ggtitle("") & NoLegend() & NoAxes()

NA
NA

Fig. 3g

dpt psuedotime in FA space

FeaturePlot(TBO_seurat, features = c("dpt_pseudotime"), pt.size=0.01,
            reduction='fa',) + scale_color_viridis(option="viridis",direction=-1)& NoAxes()

Function to plot clone dynamics colored by DPT

plotCloneDynDpt <- function(clone_pattern, seurat_obj, title = "") {
    # Subset cells of interest
    ss <- subset(seurat_obj, idents = clone_pattern)
    highlight_cells <- colnames(ss)

    # Extract embeddings and metadata
    fa_coords <- Embeddings(seurat_obj, "fa")
    dpt_vals <- seurat_obj@meta.data$dpt_pseudotime
    names(dpt_vals) <- rownames(seurat_obj@meta.data)

    dat <- as.data.frame(fa_coords)
    dat$highlight <- ifelse(rownames(dat) %in% highlight_cells, "yes", "no")
    dat$pseudotime <- dpt_vals[rownames(dat)]
    colnames(dat)[1:2] <- c("FA1", "FA2")  # name columns for clarity

    library(ggplot2)
    ggplot(dat, aes(x = FA1, y = FA2)) +
        # Background cells
        geom_point(data = subset(dat, highlight == "no"), 
                   color = "grey85", size = 0.1) +
        # Highlighted cells colored by pseudotime
        geom_point(data = subset(dat, highlight == "yes"), 
                   aes(color = pseudotime), size = 0.25) +
        scale_color_viridis_c(option = "turbo", direction = -1) +
        ggtitle(title) +
        theme_void() +
        theme(
            legend.position = "none",
            panel.border = element_blank(),      # removes any panel border
            plot.background = element_rect(fill = "transparent", color = NA),
            panel.background = element_rect(fill = "transparent", color = NA),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank()
        )
}

Fig 3h

Idents(clones) <- 'Clone_Dynamics'

my_pats <- c('Pattern_1','Pattern_2','Pattern_5','Unknown_1','Pattern_3','Pattern_4')

plots <- lapply(my_pats, function(pat) {
    tit <- gsub(pat, '_', ' ')
    plotCloneDynDpt(pat, clones, title = tit)
})

cowplot::plot_grid(plotlist = plots, ncol=6)